
## ensure you have set your working directory as appropriate (see rsr_ehsi1_loadprocess.R)

# if not already run: 
source("rsr_ehsi_2mainAnalysis.R")

#### install some additional packages:
library("car")
library("dplyr")
library("ggplot2")

##################### Figure A1: distribution of subjective job loss probabilities
par(mar = c( 5.1, 4.1, 4.1, 4.1))
barplot(summary(as.factor(ehsi$losejob))[1:11]/length(which(is.na(ehsi$losejob)==F)), names.arg = seq(0, 100, length = 11), xlab = "Subjective probability of job loss")
### Print command used for the figures in the paper (not run):
# dev.print(device = jpeg, width = 500, file = "~/Documents/research/inProgress/philosophical/risk/subj_risk_distrib.jpg")

##################### Table A5: Full results of logistic regression models of self-respect. The figures in the main text are created from models 1, 3, and 5, respectively

### Models 1, 3, and 5 in the appendix table are the models ending 2.0, run in the main paper analysis file. 
### Models 2, 4, and 6 include controls for education and marital status, and model 2 also controls for occupation; these are suffixed 2.1 and run here:

lj.sr2.1 = glm(sr.binary2 ~ losejob + topqual2 + eqinc.k + SOC2010B + selfemp + female + marital + age + I(age^2) , data = ehsi[which(ehsi$lookwk == 0 & ehsi$retired == 0 & ehsi$out_dis == 0 & ehsi$out_home == 0),], family = binomial(link = "logit"))

ue.sr2.1 = glm(sr.binary2 ~ occ_gen_uer + topqual2 + eqinc.k + selfemp + female + marital + age + I(age^2) , data = ehsi[which(ehsi$lookwk == 0 & ehsi$retired == 0 & ehsi$out_dis == 0 & ehsi$out_home == 0),], family = binomial(link = "logit"))

sk.sr2.1 = glm(sr.binary2 ~ sksp.rel + topqual2 + eqinc.k + selfemp + female + marital + age + I(age^2) , data = ehsi[which(ehsi$lookwk == 0 & ehsi$retired == 0 & ehsi$out_dis == 0 & ehsi$out_home == 0),], family = binomial(link = "logit"))

#### Table A5:
stargazer(lj.sr2.0, lj.sr2.1, ue.sr2.0, ue.sr2.1, sk.sr2.0, sk.sr2.1, type = "text", style = "ajps", omit = c("Constant", "topqual", "SOC2010", "marital", "topqual"), digits = 2, keep.stat = c("n", "aic", "ll"), order = c("losejob", "occ_gen_uer", "sksp", "eqinc", "selfemp", "female", "age"), covariate.labels = c("Subjective likelihood of job loss", "Occupational unemployment", "Skill specificity", "Equivalised income (000s)", "Self-employed", "Female", "Age", "Age 2"), add.lines = list(c("Controls: education, marital status", "No", "Yes", "No", "Yes", "No", "Yes"), c("Controls: occupation", "No", "Yes", "No", "No", "No", "No")))

##################### Table A6: Full results of logistic regression models of self-respect using a less demanding cutoff for the measurement of self-respect
lj.sr2.0b = glm(sr.binary ~ losejob + eqinc.k + selfemp + female  + age + I(age^2) , data = ehsi[which(ehsi$lookwk == 0 & ehsi$retired == 0 & ehsi$out_dis == 0 & ehsi$out_home == 0),],  family = binomial(link = "logit"))

lj.sr2.1b = glm(sr.binary ~ losejob + topqual2 + eqinc.k + SOC2010B + selfemp + female + marital + topqual2 + age + I(age^2) , data = ehsi[which(ehsi$lookwk == 0 & ehsi$retired == 0 & ehsi$out_dis == 0 & ehsi$out_home == 0),], family = binomial(link = "logit"))

ue.sr2.0b = glm(sr.binary ~ occ_gen_uer + eqinc.k + selfemp + female  + age + I(age^2) , data = ehsi[which(ehsi$lookwk == 0 & ehsi$retired == 0 & ehsi$out_dis == 0 & ehsi$out_home == 0),],  family = binomial(link = "logit"))

ue.sr2.1b = glm(sr.binary ~ occ_gen_uer + topqual2 + eqinc.k + selfemp + female + marital + topqual2 + age + I(age^2) , data = ehsi[which(ehsi$lookwk == 0 & ehsi$retired == 0 & ehsi$out_dis == 0 & ehsi$out_home == 0),], family = binomial(link = "logit"))

sk.sr2.0b = glm(sr.binary ~ sksp.rel + eqinc.k + selfemp + female  + age + I(age^2) , data = ehsi[which(ehsi$lookwk == 0 & ehsi$retired == 0 & ehsi$out_dis == 0 & ehsi$out_home == 0),],  family = binomial(link = "logit"))

sk.sr2.1b = glm(sr.binary ~ sksp.rel + topqual2 + eqinc.k + selfemp + female + marital + topqual2 + age + I(age^2) , data = ehsi[which(ehsi$lookwk == 0 & ehsi$retired == 0 & ehsi$out_dis == 0 & ehsi$out_home == 0),], family = binomial(link = "logit"))

stargazer(lj.sr2.0b, lj.sr2.1b, ue.sr2.0b, ue.sr2.1b, sk.sr2.0b, sk.sr2.1b, type = "text", style = "ajps", omit = c("Constant", "topqual", "SOC2010", "marital", "topqual2"), digits = 2, keep.stat = c("n", "aic", "ll"), order = c("losejob", "occ_gen_uer", "sksp", "eqinc", "selfemp", "female", "age"), covariate.labels = c("Subjective likelihood of job loss", "Occupational unemployment", "Skill specificity", "Equivalised income (000s)", "Self-employed", "Female", "Age", "Age 2"), add.lines = list(c("Controls: education, marital status", "No", "Yes", "No", "Yes", "No", "Yes"), c("Controls: occupation", "No", "Yes", "No", "No", "No", "No")))


##################### Table A7: : Linear models of a continuous measure of self-respect
lj.sr2.0c = lm(sr.binary ~ losejob + eqinc.k + selfemp + female  + age + I(age^2) , data = ehsi[which(ehsi$lookwk == 0 & ehsi$retired == 0 & ehsi$out_dis == 0 & ehsi$out_home == 0),])

lj.sr2.1c = lm(sr.binary ~ losejob + topqual2 + eqinc.k + SOC2010B + selfemp + female + marital + topqual2 + age + I(age^2) , data = ehsi[which(ehsi$lookwk == 0 & ehsi$retired == 0 & ehsi$out_dis == 0 & ehsi$out_home == 0),])

ue.sr2.0c = lm(sr.binary ~ occ_gen_uer + eqinc.k + selfemp + female  + age + I(age^2) , data = ehsi[which(ehsi$lookwk == 0 & ehsi$retired == 0 & ehsi$out_dis == 0 & ehsi$out_home == 0),])

ue.sr2.1c = lm(sr.binary ~ occ_gen_uer + topqual2 + eqinc.k + selfemp + female + marital + topqual2 + age + I(age^2) , data = ehsi[which(ehsi$lookwk == 0 & ehsi$retired == 0 & ehsi$out_dis == 0 & ehsi$out_home == 0),])

sk.sr2.0c = lm(sr.binary ~ sksp.rel + eqinc.k + selfemp + female  + age + I(age^2) , data = ehsi[which(ehsi$lookwk == 0 & ehsi$retired == 0 & ehsi$out_dis == 0 & ehsi$out_home == 0),])

sk.sr2.1c = lm(sr.binary ~ sksp.rel + topqual2 + eqinc.k + selfemp + female + marital + topqual2 + age + I(age^2) , data = ehsi[which(ehsi$lookwk == 0 & ehsi$retired == 0 & ehsi$out_dis == 0 & ehsi$out_home == 0),])

stargazer(lj.sr2.0c, lj.sr2.1c, ue.sr2.0c, ue.sr2.1c, sk.sr2.0c, sk.sr2.1c, type = "text", style = "ajps", omit = c("Constant", "topqual", "SOC2010", "marital", "topqual2"), digits = 2, keep.stat = c("n", "rsq", "adj.rsq"), order = c("losejob", "occ_gen_uer", "sksp", "eqinc", "selfemp", "female", "age"), covariate.labels = c("Subjective likelihood of job loss", "Occupational unemployment", "Skill specificity", "Equivalised income (000s)", "Self-employed", "Female", "Age", "Age 2"), add.lines = list(c("Controls: education, marital status", "No", "Yes", "No", "Yes", "No", "Yes"), c("Controls: occupation", "No", "Yes", "No", "No", "No", "No")))


##################### Table A8: The relationship between risk and the `first prong’: self-worth only, or the three items excluding “feeling confident”.

## Table includes original model, lj.sr2.0 as model 1

## model using only "worth"
lj.worth2.0 = glm(worth.binary2 ~ losejob + eqinc.k + selfemp + female  + age + I(age^2) , data = ehsi[which(ehsi$lookwk == 0 & ehsi$retired == 0 & ehsi$out_dis == 0 & ehsi$out_home == 0),],  family = binomial(link = "logit"))

## model using only self-respect but excluding "confidence"
lj.exconf2.0 = glm(exconf.binary2 ~ losejob + eqinc.k + selfemp + female  + age + I(age^2) , data = ehsi[which(ehsi$lookwk == 0 & ehsi$retired == 0 & ehsi$out_dis == 0 & ehsi$out_home == 0),],  family = binomial(link = "logit"))

stargazer(lj.sr2.0, lj.worth2.0, lj.exconf2.0, type = "text", style = "ajps", omit = c("eqinc", "selfemp", "female", "age", "Constant"), covariate.labels = "Subjective likelihood of job loss", notes.append = T, notes = c("All models include controls for income,  self-employment, gender, age, age squared"), digits = 2, digits.extra = 0)


##################### Table A9: Skill specificity and unemployment risk modelled simultaneously

ue.sk.sr2.1 = glm(sr.binary2 ~ occ_gen_uer + sksp.rel + eqinc.k + selfemp + female + age + I(age^2) , data = ehsi[which(ehsi$lookwk == 0 & ehsi$retired == 0 & ehsi$out_dis == 0 & ehsi$out_home == 0),], family = binomial(link = "logit"))

ue.sk.sr2.1b = glm(sr.binary ~ occ_gen_uer + sksp.rel + eqinc.k + selfemp + female + age + I(age^2) , data = ehsi[which(ehsi$lookwk == 0 & ehsi$retired == 0 & ehsi$out_dis == 0 & ehsi$out_home == 0),], family = binomial(link = "logit"))

stargazer(ue.sk.sr2.1, ue.sk.sr2.1b, type = "text", style = "ajps", omit = c("eqinc", "selfemp", "female", "age", "Constant"), covariate.labels = c("Occupational unemployment", "Skill specificity"), notes.append = T, notes = c("All models include controls for income, self-employment, gender, age, age squared"), digits = 2, digits.extra = 0)

##################### Figure A2: Marginal effects of subjective probability of job loss, modeled as a factor. 

f.lj = glm(sr.binary2 ~ losejob.f + + topqual2 + 
     eqinc.k +# lookwk + retired + out_dis + out_home + 
     selfemp + 
     as.factor(marital) + female + age + I(age^2),  
 	 family = binomial(link = "logit"), 
     data = ehsi[which(ehsi$lookwk == 0 & ehsi$retired == 0 & ehsi$out_dis == 0 & ehsi$out_home == 0),])

sue.levs = c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1)
sue.levnames = paste("losejob.f", sue.levs, sep = "")

cimat.nlf = matrix(NA, nrow = 10, ncol = 2)
for(i in 1:10){
	cimat.nlf[i,] = confint(f.lj, sue.levnames[i])
}

par(mar = c( 5.1, 4.1, 4.1, 4.1))
b=barplot(table(ehsi$losejob.f), col = rgb(0, 0.6, 0.6, 0.2), border = rgb(0, 0.6, 0.6, 0.2), xlab = "", ylab = "", yaxt = "n", main = "", xaxt = "n", ylim = c(0, 1200))
axis(4, col = rgb(0, 0.6, 0.6, 0.3), col.ticks = rgb(0, 0.6, 0.6, 0.3), col.axis = rgb(0, 0.6, 0.6, 0.3), col.lab = rgb(0, 0.6, 0.6, 0.3),at = c(200, 400, 600, 800, 1000))
mtext("Number of respondents", side = 4, line = 3, col = rgb(0, 0.6, 0.6, 0.3))
par(new = T)
plot(c(0,sue.levs), c(0,coefficients(f.lj)[2:11]), ylim = c(-2,2), xlab = "Subjective Probability of Job Loss", ylab = "Marginal effect of change of risk on linear predictor", xlim = c(-0.04, 1.04))
lines(c(0,sue.levs), c(0,coefficients(f.lj)[2:11]), col = rgb(0, 0.6, 0.6, 0.7), lwd = 2)
segments(x0 = sue.levs, x1 = sue.levs, y0 = cimat.nlf[,1], y1 = cimat.nlf[,2])
abline(h= 0, col = "grey60")

#### Print command used for the figures in the paper (not run):
# dev.print(device = jpeg, width = 500, file = "factor_lj_mfx_withHist.jpg")

##################### Table A10: Non-linear models, subjective job loss

# linear (original)
l.lj = glm(sr.binary2 ~  losejob + topqual2 + 
     eqinc.k +# lookwk + retired + out_dis + out_home + 
     selfemp + 
     as.factor(marital) + female + age + I(age^2),  
	 family = binomial(link = "logit"), 
     data = ehsi[which(ehsi$lookwk == 0 & ehsi$retired == 0 & ehsi$out_dis == 0 & ehsi$out_home == 0),])

# quadratic
q.lj = glm(sr.binary2 ~ losejob + I(losejob^2) + topqual2 + 
     eqinc.k +# lookwk + retired + out_dis + out_home + 
     selfemp + 
     as.factor(marital) + female + age + I(age^2),  
 	 family = binomial(link = "logit"), 
     data = ehsi[which(ehsi$lookwk == 0 & ehsi$retired == 0 & ehsi$out_dis == 0 & ehsi$out_home == 0),])

# step at 0.3
ehsi$lj0.3 = as.numeric(ehsi$losejob >= 0.3)
l.lj.f3 = glm(sr.binary2 ~ lj0.3 + topqual2 + 
     eqinc.k +# lookwk + retired + out_dis + out_home + 
     selfemp + 
     as.factor(marital) + female + age + I(age^2),  
 	 family = binomial(link = "logit"), 
     data = ehsi[which(ehsi$lookwk == 0 & ehsi$retired == 0 & ehsi$out_dis == 0 & ehsi$out_home == 0),])


stargazer(l.lj, q.lj, l.lj.f3, f.lj, omit = c("Constant", "female", "topqual", "eqinc", "age", "marital", "selfemp"), type = "text", dep.var.labels = "Binary Self-Respect", 
#dep.var.caption = c("Linear", "Quadratic", "Step at 0.3", "Factor"), 
 style = "ajps", covariate.labels = c("Subjective probability of job loss", "Subjective probability squared", "Subjective probability > 0.3", "Subjective probability = 0.1", "Subjective probability = 0.2", "Subjective probability = 0.3", "Subjective probability = 0.4", "Subjective probability = 0.5", "Subjective probability = 0.6", "Subjective probability = 0.7", "Subjective probability = 0.8", "Subjective probability = 0.9", "Subjective probability = 1"))


##################### Table A11: Non-linear models, objective risk measures
 l.ur2 = glm(sr.binary2 ~ occ_gen_uer + topqual2 + 
     eqinc.k +# lookwk + retired + out_dis + out_home + 
     selfemp + 
     as.factor(marital) + female + age + I(age^2),  
 	 family = binomial(link = "logit"), 
     data = ehsi[which(ehsi$lookwk == 0 & ehsi$retired == 0 & ehsi$out_dis == 0 & ehsi$out_home == 0),])

 q.ur2 = glm(sr.binary2 ~ occ_gen_uer + I(occ_gen_uer^2) + topqual2 + 
     eqinc.k +# lookwk + retired + out_dis + out_home + 
     selfemp + 
     as.factor(marital) + female + age + I(age^2),  
 	 family = binomial(link = "logit"), 
     data = ehsi[which(ehsi$lookwk == 0 & ehsi$retired == 0 & ehsi$out_dis == 0 & ehsi$out_home == 0),])
 
 
 l.sksp2 = glm(sr.binary2 ~ sksp.rel + topqual2 + 
     eqinc.k +# lookwk + retired + out_dis + out_home + 
     selfemp + 
     as.factor(marital) + female + age + I(age^2),  
 	 family = binomial(link = "logit"), 
     data = ehsi[which(ehsi$lookwk == 0 & ehsi$retired == 0 & ehsi$out_dis == 0 & ehsi$out_home == 0),])

 q.sksp2 = glm(sr.binary2 ~ sksp.rel + I(sksp.rel^2) + topqual2 + 
     eqinc.k +# lookwk + retired + out_dis + out_home + 
     selfemp + 
     as.factor(marital) + female + age + I(age^2),  
 	 family = binomial(link = "logit"), 
     data = ehsi[which(ehsi$lookwk == 0 & ehsi$retired == 0 & ehsi$out_dis == 0 & ehsi$out_home == 0),])

 stargazer(l.ur2, q.ur2, l.sksp2, q.sksp2, omit = c("Constant", "female", "topqual", "eqinc", "age", "marital", "selfemp"), type = "text", dep.var.labels = "Binary Self-Respect", style = "ajps", covariate.labels = c("Occupational unemployment", "Occupational unemployment squared", "Skill specificity", "Skill specificity squared"), order = c("occ_", "sksp"), notes.append = T, notes = c("All models include controls for education, income, self-employment, marital status, education, gender, age, and age squared."), title = "Table A11: Non-linear effects of risk: objective measures. Standard errors in parentheses.")
  
##################### Figure A3: The effect of objective unemployment risk on self-respect. Marginal effects estimated from model 2 (table A10)

 z0.ur <- seq(min(ehsi$occ_gen_uer, na.rm=T), max(ehsi$occ_gen_uer, na.rm=T), length.out = 1000)
## marginal effects of changes in risk: 
beta.hat.ur = coef(q.ur2)
dy.dx.ur <- beta.hat.ur["occ_gen_uer"]*z0.ur + beta.hat.ur["I(occ_gen_uer^2)"]*z0.ur^2

beta_ur <- coef(q.ur2)["occ_gen_uer"]
beta_ur_2 <- coef(q.ur2)["I(occ_gen_uer^2)"]

cov.ur = vcov(q.ur2)

se.dy.dx.ur <- sqrt(cov.ur["occ_gen_uer", "occ_gen_uer"] + z0.ur^2*cov.ur["I(occ_gen_uer^2)", "I(occ_gen_uer^2)"] + 2*z0.ur*cov.ur["occ_gen_uer", "I(occ_gen_uer^2)"])
 # calculate upper and lower bounds of a 95% CI
upr.ur <- dy.dx.ur + 1.96*se.dy.dx.ur
lwr.ur <- dy.dx.ur - 1.96*se.dy.dx.ur

 plot(z0.ur, dy.dx.ur, type = "l", lwd = 2, xlab = "Occupation-Gender Unemployment Risk", ylab = "Marginal effect of change of risk on linear predictor", ylim = c(-1.5, 0.05), lty = 2)
 
 polygon(c(z0.ur, rev(z0.ur)), c(upr.ur, rev(lwr.ur)), border = rgb(0, 0, 0, 0.1), col = rgb(0, 0, 0, 0.1))
 
### Print command used for the figures in the paper (not run):
# dev.print(device = jpeg, width = 500, file = "ehsi_uer_quadratic.jpg")

##################### Figure A4: The effect of objective skill-specificty risk on self-respect. Marginal effects estimated from model 4 (table A10)

z0.sksp <- seq(min(ehsi$sksp.rel, na.rm=T), max(ehsi$sksp.rel, na.rm=T), length.out = 1000)
beta.hat.sksp = coef(q.sksp2)
dy.dx.sksp <- beta.hat.sksp["sksp.rel"]*z0.sksp + beta.hat.sksp["I(sksp.rel^2)"]*z0.sksp^2

beta_sksp <- coef(q.sksp2)["sksp.rel"]
beta_sksp_2 <- coef(q.sksp2)["I(sksp.rel^2)"]

cov.sksp = vcov(q.sksp2)

se.dy.dx.sksp <- sqrt(cov.sksp["sksp.rel", "sksp.rel"] + z0.sksp^2*cov.sksp["I(sksp.rel^2)", "I(sksp.rel^2)"] + 2*z0.sksp*cov.sksp["sksp.rel", "I(sksp.rel^2)"])
 # calculate upper and lower bounds of a 95% CI
upr.sksp <- dy.dx.sksp + 1.96*se.dy.dx.sksp
lwr.sksp <- dy.dx.sksp - 1.96*se.dy.dx.sksp

 plot(z0.sksp, dy.dx.sksp, type = "l", lwd = 2, xlab = "Skill Specificity", ylab = "Marginal effect of change of risk on linear predictor", ylim = c(-0.5, 0.1), lty = 2)
# abline(h=0, col = "grey60")
 polygon(c(z0.sksp, rev(z0.sksp)), c(upr.sksp, rev(lwr.sksp)), border = rgb(0, 0, 0, 0.1), col = rgb(0, 0, 0, 0.1))
 
 ### Print command used for the figures in the paper (not run):
# dev.print(device = jpeg, width = 500, file = "ehsi_sksp_quadratic.jpg")

##################### Figure A5: the effect of risk on self-respect at different income levels. Estimates generated from model 1 (table A12). 

lj.incbufsr2.0 = glm(sr.binary2 ~ losejob*eqinc.k + selfemp + female  + age + I(age^2) , data = ehsi[which(ehsi$lookwk == 0 & ehsi$retired == 0 & ehsi$out_dis == 0 & ehsi$out_home == 0),],  family = binomial(link = "logit"))

ld <- seq(0, 1, 0.025)

pry14 = predict(lj.incbufsr2.0, 
	data.frame(losejob = ld,
   	eqinc.k = rep(14, length(ld)),
   	selfemp = rep(0, length(ld)), 
   	female = rep(0, length(ld)), 
   	age = rep(45, length(ld))
   	), type = "response", se.fit = T)

pry41 = predict(lj.incbufsr2.0, 
	data.frame(losejob = ld,
   	eqinc.k = rep(41, length(ld)),
   	selfemp = rep(0, length(ld)), 
   	female = rep(0, length(ld)), 
   	age = rep(45, length(ld))
   	), type = "response", se.fit = T)

lower14 = pry14$fit - 1.96*pry14$se.fit
upper14 = pry14$fit + 1.96*pry14$se.fit

lower41 = pry41$fit - 1.96*pry41$se.fit
upper41 = pry41$fit + 1.96*pry41$se.fit


plot(c(0,1), c(0,0.7), type = "n", xlab = "Subjective Probability of Job Loss", ylab = "Predicted Probability of Reporting Self-Respect")

# se color in cplot: grDevices::gray(0.5, 0.5)

polygon(c(ld, rev(ld)), c(lower14, rev(upper14)), col = grDevices::gray(0.5, 0.5), border = grDevices::gray(0.5, 0.5))

polygon(c(ld, rev(ld)), c(lower41, rev(upper41)), col = grDevices::gray(0.5, 0.5), border = grDevices::gray(0.5, 0.5))


lines(ld, pry14$fit, , lty = 2)
lines(ld, pry41$fit)

legend(0, 0.2, bty = "n", lty = c(1, 2),  legend = c("High income (75th percentile)", "Low income (25th percentile)"))
legend(0.03, 0.09, bty = "n", fill = grDevices::gray(0.5, 0.5), legend = "95% confidence intervals", border = grDevices::gray(0.5, 0.5))

### Print command used for the figures in the paper (not run):
# dev.print(device = jpeg, width = 500, file = "subjR_incomebuffer.jpg")


##################### Figure A6: the effect of risk on self-respect with and without partner. Estimates generated from model 2 (table A12). 
ehsi$part = as.numeric(ehsi$marital == 1)
ehsi$part = as.factor(ehsi$part)

lj.pbufsr2.0 = glm(sr.binary2 ~ losejob*part + selfemp + female + age + I(age^2) , data = ehsi[which(ehsi$lookwk == 0 & ehsi$retired == 0 & ehsi$out_dis == 0 & ehsi$out_home == 0),],  family = binomial(link = "logit"))

prb0 = predict(lj.pbufsr2.0, 
	data.frame(losejob = ld,
   	part = rep("0", length(ld)),
   	#marital = rep("0", length(ld)),
   	selfemp = rep(0, length(ld)), 
   	female = rep(0, length(ld)), 
   	age = rep(45, length(ld)), 
	eqinc.k = rep(mean(ehsi$eqinc.k, na.rm = T), length(ld))
   	), type = "response", se.fit = T)

prb1 = predict(lj.pbufsr2.0, 
	data.frame(losejob = ld,
   	part = rep("1", length(ld)),
   #	marital = rep("1", length(ld)),
   	selfemp = rep(0, length(ld)), 
   	female = rep(0, length(ld)), 
   	age = rep(45, length(ld)), 
   	eqinc.k = rep(mean(ehsi$eqinc.k, na.rm = T), length(ld))
   	), type = "response", se.fit = T)

lower1 = prb1$fit - 1.96*prb1$se.fit
upper1 = prb1$fit + 1.96*prb1$se.fit

lower0 = prb0$fit - 1.96*prb0$se.fit
upper0 = prb0$fit + 1.96*prb0$se.fit


plot(c(0,1), c(0,0.7), type = "n", xlab = "Subjective Probability of Job Loss", ylab = "Predicted Probability of Reporting Self-Respect")

# se color in cplot: grDevices::gray(0.5, 0.5)

polygon(c(ld, rev(ld)), c(lower1, rev(upper1)), col = grDevices::gray(0.5, 0.5), border = grDevices::gray(0.5, 0.5))

polygon(c(ld, rev(ld)), c(lower0, rev(upper0)), col = grDevices::gray(0.5, 0.5), border = grDevices::gray(0.5, 0.5))


lines(ld, prb0$fit, , lty = 2)
lines(ld, prb1$fit)

legend(0, 0.2, bty = "n", lty = c(1, 2),  legend = c("Married/cohabiting", "Single/divorced/widowed"))
legend(0.03, 0.09, bty = "n", fill = grDevices::gray(0.5, 0.5), legend = "95% confidence intervals", border = grDevices::gray(0.5, 0.5))

### Print command used for the figures in the paper (not run):
# dev.print(device = jpeg, width = 500, file = "subjR_partnerbuffer.jpg")

##################### Table A12: The effect of risk on self-respect for different kinds of respondent
stargazer(lj.incbufsr2.0, lj.pbufsr2.0, type = "text", 
	omit = c("self", "female", "age", "Constant"), order = c(":", "losejob","eqinc", "part1"), covariate.labels = c("Job loss probability x income", "Job loss probability x adult partner", "Subjective probability of job loss", "Income (000s)","Adult partner in household"),
	notes.append = T, notes = c("Models include controls for gender, age, and age squared."),
	title = "Table A12: The effect of risk on self-respect for different types of respondent. Standard errors in parentheses.",
	digits = 2, digits.extra = 0,  style = "ajps")


##################### Table A13: The effect of risk exposure on self-respect, for varying levels of liking risk
lj.sr2b.risk = glm(sr.binary2 ~ losejob*riskseek + eqinc.k + selfemp + female  + marital + topqual2 + age + I(age^2) , data = ehsi[which(ehsi$lookwk == 0 & ehsi$retired == 0 & ehsi$out_dis == 0 & ehsi$out_home == 0),],  family = binomial(link = "logit"))

ue.sr2b.risk = glm(sr.binary2 ~ occ_gen_uer*riskseek + eqinc.k + selfemp + female  + marital + topqual2 + age + I(age^2) , data = ehsi[which(ehsi$lookwk == 0 & ehsi$retired == 0 & ehsi$out_dis == 0 & ehsi$out_home == 0),],  family = binomial(link = "logit"))

sk.sr2b.risk = glm(sr.binary2 ~ sksp.rel*riskseek + eqinc.k + selfemp + female  + marital + topqual2 + age + I(age^2) , data = ehsi[which(ehsi$lookwk == 0 & ehsi$retired == 0 & ehsi$out_dis == 0 & ehsi$out_home == 0),],  family = binomial(link = "logit"))

stargazer(lj.sr2b.risk, ue.sr2b.risk, sk.sr2b.risk, type = "text", style = "ajps", omit = c("Constant", "topqual", "SOC2010", "marital", "topqual"), digits = 2, digits.extra = 0, keep.stat = c("n", "aic", "ll"), order = c("riskseek",":", "losejob", "occ_gen_uer", "sksp","eqinc", "selfemp", "female", "age"), dep.var.labels = "Binary self-respect", covariate.labels = c("Liking risk", "Likelihood of job loss x liking risk", "Occupational unemployment x liking risk", "Skill specificity x liking risk", "Subjective probability of job loss", "Occupational unemployment rate", "Skill specificity", "Equivalised income (000s)", "Self-employed", "Female", "Age", "Age 2"), notes.append = T, notes = c("Models include controls for education, income, gender, marital status, age, and age squared."))

